#include <iostream>
#include <string>
#include <Magick++.h>

#include "../common/tools.h"
#include "../common/performance.h"
#include "../common/imgtools.h"

#include "LBM.h"
#include "visualizers.h"
#include "particles.h"


using std::cout;
using std::endl;
using std::string;

using namespace Magick;

const LBfloat FlowVelX = 0.1;
const LBfloat FlowVelY = 0.0;
const LBfloat FlowDensity = 4.5;
const LBfloat FlowViscosity = 1.0 / 50.0;

const int GridWidth = 512;
const int GridHeight = 300;
const int ObstacleRadius = 20;

const int MaxIteration = 20000;
const int DumpInterval = 50;

const LBfloat ParticleDensity = 0.03;

LatticeBoltzmannD2Q7 ca(GridWidth, GridHeight);
ParticleSystem particles(ca);


void dumpCA(int iterNum)
{
  Image img = renderCA(ca, VelHSLPainter(8.0));
  
  int width = ca.getWidth(), height = ca.getHeight();
  PixelPacket *pixel_cache = img.getPixels(0, 0, width, height);

  ParticleSystem::const_iterator iter = particles.begin();
  while ( iter != particles.end() )
  {
    int xi = (int) iter->x;
    int yi = (int) iter->y;
    if (xi >=0 && xi < width && yi >= 0 && yi < height)
    {
      PixelPacket & p = *(pixel_cache + yi * width + xi);
      ColorRGB c(p);
      p = ColorRGB(1.0 - c.red(), 1.0 - c.green(), 1.0 - c.blue());
    }
    ++iter;
  }

  string iterStr = format("{0}") % iterNum;
  string zeros = string(4 - iterStr.length(), '0');
  iterStr = zeros.append(iterStr);
  img.write( format("dump_{0}.jpg") % iterStr );
}

void fillTypeEllips(LatticeBoltzmannD2Q7 & ca, int x, int y, int w, int h, CellType type)
{
  for (int cy = 0; cy != h; ++cy)
    for (int cx = 0; cx != w; ++cx)
    {
      float dx = (float) 2.0f * cx / (w) - 1.0f;
      float dy = (float) 2.0f * cy / (h) - 1.0f;
      if (dx*dx + dy*dy <= 1.0f)
        ca.setCellType( x+cx, y+cy, type);
    }
}

void initCA()
{
  ca.setViscosity(FlowViscosity);
  D2Q7Cell eq;
  ca.calcEquilibrium(FlowDensity, FlowVelX, FlowVelY, eq);
  ca.setConstBound(CELL_EMPTY, eq);

  for (int y = 0; y != GridHeight; ++y)
    for (int x = 0; x != GridWidth; ++x)
      ca.setCell(x, y, eq);

  int w = GridWidth;
  int h = GridHeight;
  int rw = ObstacleRadius;
  int rh = (int) (ObstacleRadius / CellHeight);
  fillTypeEllips(ca, w / 5 - rw, h / 2 - rh, 2*rw, 2*rh, CELL_WALL);
}

int main()
{
  initCA();

  Image velMap = renderVelocityMap(100);
  velMap.write("vel_map.jpg");

  size_t particleCount = (int) GridWidth * GridHeight * ParticleDensity;
  for (int i = 0; i != particleCount; ++i)
    particles.addParticle(randf() * GridWidth, randf() * GridHeight);
  particleCount = particles.getCount();

  int dumpCount = 0;
  PerformanceTimer timer;
  timer.start();
  do {
    timer.start();
    ca.iterate();
    timer.stop();
    particles.iterate();
    while (particles.getCount() < particleCount)
      particles.addParticle(0, randf() * GridHeight);

    if ( ca.getIteration() % DumpInterval == 0)
    {
      cout << "IPS: " << 1.0 / timer.getMeanTime() << endl;
      timer.reset();

      cout << "saving iteration " << ca.getIteration() << " ... ";
      dumpCA(++dumpCount);
      cout << "ok" << endl;
    }
  } while ( ca.getIteration() != MaxIteration );

  return 0;
}
